Computational Investigation of Interactions between Carbon Nitride Dots and Doxorubicin

The study of carbon dots is one of the frontiers of materials science due to their great structural and chemical complexity. These issues have slowed down the production of solid models that are able to describe the chemical and physical features of carbon dots. Recently, several studies have started to resolve this challenge by producing the first structural-based interpretation of several kinds of carbon dots, such as graphene and polymeric ones. Furthermore, carbon nitride dot models established their structures as being formed by heptazine and oxidized graphene layers. These advancements allowed us to study their interaction with key bioactive molecules, producing the first computational studies on this matter. In this work, we modelled the structures of carbon nitride dots and their interaction with an anticancer molecule (Doxorubicin) using semi-empirical methods, evaluating both geometrical and energetic parameters.


Introduction
Since their discovery in 2004 [1], carbon nano dots (CDs) have attracted great interest due to their unique set of properties [2]. CDs are characterized by a high fluorescence emission [3] due to their quantum confinement [4], water solubility [5] and high biocompatibility [6]. Their chemical features have spread their use across several fields with remarkable achievements as fluorescent probes in analytical procedures [7][8][9] and biomedical treatments [10,11] as active materials for both photocatalysis [12,13] and electrocatalysis [14,15]. The application of CDs has raised several issues related to what is and what is not a carbon dot together with serious concern about the development of a rational classification methodology [2].
The classification is of utmost relevance to properly approaching the field of CDs [16]. The most widespread classification approach is based on structural features, and it encompasses three large families of CDs [17]: (i) graphene quantum dots (GQDs), (ii) carbon nitride carbon dots (CNDs) and (iii) polymeric carbon dots (PCDs).
GQDs represent the first kind of CDs discovered and they present the simplest structures among all CDs. The chemical structure of this family is directly related to small graphene oxide clusters, and the distribution of chemical functionalities is explained simply by the Lerf-Klinowsky model [18]. Accordingly, GQDs contain few layers of oxidized graphene, with carbonyls and carboxylic residues being concentrated in the edges, producing an oxygen-rich shell that accounts for their water solubility. The oxidized graphene layers of GQDs display different morphologies and oxidation degrees [19,20] that provide different chemical platforms for further functionalization [21,22].
PCDs are a wide family of nanomaterials produced by the partial degradation of polymers [23] or through the condensation of organic precursors [24][25][26]. Their properties are very difficult to correlate with a unique interpretative model and each PCD requires proper characterization.
CNDs are similar to GQDs, but they present a more complex structure in which nitrogen atoms dope the graphene layers [27] or form nitrogen-rich aromatic moieties, as found by Mintz et al. [28]. This represents a very interesting possibility for synthetic chemistry, as it can tune the reactivity and properties [29], particularly in a watery medium, as mentioned by Wiśniewski [30]. As a consequence of their functional tuneability, CNDs can act as a platform for several specific tasks, such as chemotherapy for cancer treatment, by being conjugated with active molecules [31]. Among them, Doxorubicin is one of the most studied, due to its remarkable activity, biological stability and simple chemical structure [32,33]. CNDs were successfully combined with Doxorubicin for the production of conjugates (D-CNDs) for the treatment of several cancers using theragnostic approaches [34][35][36]. Even if the interest in the rational design of D-CNDs is great, the literature lacks computational studies devoted to quantifying the interactions between Doxorubicin and CNDs. Rashid et al. [37] investigated flutamide CNDs through DFT calculation, limiting the system to a representative portion of small CNDs. Zaboli et al. [38] went more in depth by simulating the interaction of a molecule of Doxorubicin on a single sheet of carbon nitride, taken as a representative model of CNDs. The authors were able to correctly evaluate the π−π and water interaction with both a single large layer of heptazine and the Doxorubicin-conjugated material. Similarly, Shaki et al. [39] investigated the adsorption of and interaction between several anticancer species with the external and inner surfaces of carbon nanotubes. All of these studies show a common point represented by the charge transfer from the drugs to the aromatic domain-rich carrier and assume that this is the key reason for the efficacy of the systems.
These approaches are powerful, but they require a high calculation power, and they are currently limited to small chemical systems that are not sufficiently complex to properly describe CDs. Alternatively, semi-empirical quantum chemistry has been developed to resolve the computation challenges of complex systems in order to try to overcome the main limitations of the traditional Hartree-Fock approaches, the slow speed and the low accuracy [40]. The semi-empirical methods are based on the time reduction of the calculation by parameterizing or omitting certain terms based on the experimental data sources (i.e., the ionization energies of atoms, dipole moments of molecules) [41]. Accordingly, the semi-empirical quantum chemistry approaches are considerably faster and useful for modelling large molecules. The first semi-empirical quantum chemistry method was the modified neglect of the diatomic overlap (MNDO) [42] that parameterized the one-center two-electron integrals based on the spectroscopic data for the isolated atoms, evaluating other two-electron integrals using the multipole-multipole interactions of the classical electrostatics [43]. The MNDO approach was implemented using a basis set composed of only s and p orbitals even if the d orbitals set was used to describe hypervalent sulfur atoms and transition metals. Nevertheless, MNDO is characterized by several intrinsic limitations, such as the impossibility to describe both hydrogens or to predict the heat of the formation [44]. Dewar and coworkers [45] improved the MNDO approach by developing the Austin Model X (AMX) by modifying the expression of nuclear-nuclear core repulsion, emulating van der Waals interactions. This new approach required a total reparameterization of the dipole moments and ionization potentials but allowed the description of the hydrogen bond network and the heat of the formation without properly estimating the basicity.
All of these approaches were made outdated by the development of parametric method 3 (PM3) [46], which uses a Hamiltonian operator that is very similar to the AM1 Hamiltonian, but the parameterization strategy is one of molecular properties, rather than atomic ones. A consolidating method which can be used to perform the computational modelling of large organic systems with lower computational costs is the PM3. This is a semi-empirical quantum mechanical parameterization method based on the modified neglect of the differential overlap method [47].
The PM3 approach is based on the use of a set of empirical parameters to describe the electronic structure of a molecule [48]. The key feature of the PM3 is that it offers the chance to accurately predict the geometrical conformation of large molecules, including enzymes [49] and polymers [50].
Despite these advantages, the PM3 method does not perform well enough to accurately describe the properties of the transition metal complexes [51] and the effects of electron correlation [52].
In the present work, we report a computational study focused on the evaluation of the geometrical features of D-CNDs and the interaction energies (E i ) between Doxorubicin and CNDs using a PM3-based computational approach. We ran a semi-empirical simulation of four-layer CNDs with a single molecule of Doxorubicin, evaluating different geometrical interactions with the Doxorubicin of the outer layer of the CNDs and evaluating the intercalant agent in order to provide a first solid insight into the non-covalent interactions occurring within a chemotherapy agent.

CNDs Model Structure: Definition and Modelling
The most challenging issue in the modelling of CNDs is to provide a representative species to be studied. The structural composition of CNDs is highly dependent on the synthetic method used for their production. CNDs are generally composed of a core composed of sp 2 -hybridized carbons surrounded by a surface layer that contains a mixture of sp 2 -/sp 3 -hybridized carbon and nitrogen atoms [53]. Regarding pure carbon nitride structures, the nitrogen atoms are incorporated directly into their carbon lattice, creating structural defects and altering the electronic properties of the material. Considering CNDs, the distribution of nitrogen atoms is still a matter of discussion, and several hypotheses have been suggested to explain their several spatial arrangements.
Firstly, those conducting CND structural investigations have hypothesized a pure layered structure formed by functionalized heptazine units [27]. Mintz et al. [28] moved a step forward in the clear definition of the structure of CNDs, suggesting a more complex arrangement of the layers. Based on a detailed physical-chemical investigation, the authors proved that CNDs have a graphitic core with massive functionalization on the edges and less heteroatomic inclusion in the core. As reported by Zhou et al. [54], elucidating the differences between a pure heptazine and a realistic model was crucial in order to properly evaluate the interaction between the CNDs and the protein active site. Nevertheless, the great variability of their functionalization prevented the realization of a general model compound that is able to describe all CNDs species. Accordingly, we tentatively propose a model that considers both the heptazine and functionalized graphene layers, as shown in Figure 1.
The heptazine structure (L h n ) was assembled following the method used by Zhou et al. [54] to produce CNDs composed of pure heptazine. The graphene layer L g n could be modeled using at least four geometrical models, as reported by Mandal et al. [19]. We selected type-4 as it showed the lowest energy band gap. According to the Lerf-Klinowski model of graphene oxide [18], oxygen-containing residues symmetrically tailored the L g n layer on the edges with oxygen-based functionalities. The oxygen-based functions of the CND edges are still a matter of debate, but Kirbas et al. [29] have reported how they can be tuned by varying the amounts of urea and organic acids used for their production. We selected, as oxygen-based functionalities, hydroxyl and carboxylic groups to evaluate the effects of the presence of a strong network of hydrogen bonds inside the structure. This condition could promote the dismutation of carbonyl groups under harsh synthetic environments while avoiding the presence of carbonyl residues [23]. The other key point to be defined was related to the molecular weight of L g n . We limited the L g n size to 81 carbon atoms in order to achieve a four-layer structure with two each of L h n and L g n , respectively, for an average molecular weight higher than 4000 g/mol but with an expected size of below 3 nm, in accordance with the data reported in the literature regarding the common size of CNDs [55]. The layered structure was optimized in vacuo at a temperature of 303.15 K, and the graphical results are shown in Figure 2. The heptazine structure (Lh n ) was assembled following the method used by Zhou et al. [54] to produce CNDs composed of pure heptazine. The graphene layer Lg n could be modeled using at least four geometrical models, as reported by Mandal et al. [19]. We selected type-4 as it showed the lowest energy band gap. According to the Lerf-Klinowski model of graphene oxide [18], oxygen-containing residues symmetrically tailored the Lg n layer on the edges with oxygen-based functionalities. The oxygen-based functions of the CND edges are still a matter of debate, but Kirbas et al. [29] have reported how they can be tuned by varying the amounts of urea and organic acids used for their production. We selected, as oxygen-based functionalities, hydroxyl and carboxylic groups to evaluate the effects of the presence of a strong network of hydrogen bonds inside the structure. This condition could promote the dismutation of carbonyl groups under harsh synthetic environments while avoiding the presence of carbonyl residues [23]. The other key point to be defined was related to the molecular weight of Lg n . We limited the Lg n size to 81 carbon atoms in order to achieve a four-layer structure with two each of Lh n and Lg n , respectively, for an average molecular weight higher than 4000 g/mol but with an expected size of below 3 nm, in accordance with the data reported in the literature regarding the common size of CNDs [55]. The layered structure was optimized in vacuo at a temperature of 303.15 K, and the graphical results are shown in Figure 2.  The optimized structure was characterized by a free energy of 180.4 kcal/mol and a maximum size of 2.0 × 1.3 nm, with different layer distances related to the intrinsic asymmetry of the model used. The Lg 1 -Lg 2 distance was found to be up to 0.41 nm, while the maximum distance between the Lg n and Lh n was lower: 0.39-0.34. The Lg 1 -Lg 2 interlayer distances were considerably higher than that of the neat graphite structure (0.335 nm) [56], which is in agreement with the massive functionalization that induced this defectiveness. The presence of functional groups on the edge of the Lg n induced a deformation due to the simultaneous effects of steric hindrance and the attraction/repulsion between the oxygenbased functions [57,58]. The Lh 1 and Lh 2 layer was distorted, with dihedral angles of 1.4° and 7.8° due to the interaction between the nitrogen residues and oxygen-based functionalities. Interestingly, Lh n was non-centered on the Lg n , leaving a more oxygen-rich region on one of the CNDs' extremities. This was due to the asymmetry of the Lg n , which can The optimized structure was characterized by a free energy of 180.4 kcal/mol and a maximum size of 2.0 × 1.3 nm, with different layer distances related to the intrinsic asymmetry of the model used. The L g 1 -L g 2 distance was found to be up to 0.41 nm, while the maximum distance between the L g n and L h n was lower: 0.39-0.34. The L g 1 -L g 2 interlayer distances were considerably higher than that of the neat graphite structure (0.335 nm) [56], which is in agreement with the massive functionalization that induced this defectiveness. The presence of functional groups on the edge of the L g n induced a deformation due to the simultaneous effects of steric hindrance and the attraction/repulsion between the oxygen-based functions [57,58]. The L h 1 and L h 2 layer was distorted, with dihedral angles of 1.4 • and 7.8 • due to the interaction between the nitrogen residues and oxygen-based functionalities. Interestingly, L h n was non-centered on the L g n , leaving a more oxygen-rich region on one of the CNDs' extremities. This was due to the asymmetry of the L g n , which can promote the favoring of the formation of a hydrogen bonds network on the edges, rather than π−π stacking.
As a matter of fact, this simple model was a rough approximation of CNDs, neglecting a real distribution of the oxygen and nitrogen functions between L g n and L h n , but it was worth considering it due to the simple approach proposed to introduce the grapheneoxidized core into the CNDs. We discuss the limitations and applications of our approach in the dedicated section below.

Doxo@CNDs Model Structures Modelling
The four-layer model adopted to enable us to describe the CNDs allowed several kinds of interaction with Doxorubicin, forming several D-CNDs. We investigated the structures reported in Figure 3    D-CNDs1 showed a consistent increase in terms of L h 1 -L g 1 distance compared with CNDs up to 0.55 nm with a decrement of the L h n strains and significant increase in L g n up to 15.9 • . This modification occurred due to the strong interaction of Doxorubicin with L h 1 at a distance of 0.25 nm. As it clearly emerged, the decrement of the L h n strain compacted the L g n , promoting the interactions through weak forces. When the Doxorubicin intercalated the CNDs layers, greater changes occurred in the CNDs' structure. D-CNDs2 showed an increase in all layers' distance and the strains. Interestingly, Doxorubicin was closer to L g 1 (0.36 nm) than L h 1 (0.40 nm), probably due to a better interaction with the graphene layer as a consequence of its greater spatial dimension. Similar distances were observed in D-CNDs3, where the distance between the Doxorubicin and L g n layer was the same (0.36-0.37 nm). In this case, we observed the maximum strain of L g n , reaching 63.2 • due to the formation of a task-like structure in which Doxorubicin found its place. D-CNDs4 proved that the interaction with all the CNDs layers also induced a significant structural modification, with the increase in L h 1 -L g 1 distance being up to 0.42 nm and the layers being distorted. These computational results open an interesting pathway to evaluating the microscopic analysis of D-CNDs in which structural modification could be correlated with the kind of conjugation that has occurred in the specie. Considering the values of E i , the more stable structure was D-CNDs3, in which Doxorubicin is fully contained in the CNDs structures and stabilized by the π−π interactions of L g n (E i −39.5 kcal/mol), while the D-CNDs2 and D-CNDs4 were less stable (E i 7.1 and 16.6 kcal/mol, respectively). D-CNDs1 was also stabilized by the ordering of the piling-up of Doxorubicin on L h 1 with an E i −28.5 kcal/mol. These results suggest that the π−π interactions were the main driving force for the interaction between the Doxorubicin and CNDs, rather than the hydrogen bond ones.
Doxorubicin loading to the CND also affected the HOMO-LUMO gap (∆ HL ). The calculated Doxorubicin ∆ HL was in agreement with the one reported by Lopez-Chavez et al. [59] using DFT approaches. The CNDs showed an ∆ HL of up to 1.56 eV, while the presence of Doxorubicin altered the ∆ HL by both increasing or decreased it based on its position. D-CNDs1 and D-DCNDs4 (1.38 and 1.46 eV, respectively) showed an ∆ HL lower than CNDs, while both D-CNDs2 and D-CNDs3 showed sensible increments of ∆ HL (1.71 and 1.58 eV, respectively). As reported by Bharathy et al. [60], a decrease in ∆ HL could be correlated to an increase in the reactivity of the systems, and this suggested that Doxorubicin-containing CNDs were less reactive when Doxorubicin was inserted between the layers. Particularly, the insertion between L h 1 and L g 1 seemed to be the best configuration among the ones tested due to the interaction with two layers without a massive alteration, as in the case of D-CNDs4. The insertion between L g 1 and L g 2 was not so effective, probably due to the exposure of a consistent part of Doxorubicin to the external environment. Even if ∆ HL is a powerful tool for discussing the stability of a system, it is not entirely sufficient to describe the kind of interaction occurring in the system. According to Johnson et al. [61], the evaluation of interactions could be properly conducted using only NBO analysis to evaluate the atom-atom contribution to the stabilization. Nonetheless, this approach has a very high computational cost for large systems such as the one that is very representative of those used in the interpretation of CDs. Alternatively, it is possible to evaluate the electron density and use electron maps such those reported in Figure 4. The electron density map showed a rich region between Lg 1 and Lg 2 in CNDs that was reduced in the D-CNDs, with the exception of D-CNDs2. In this case, the insertion of Doxorubicin between Lh 1 and Lg 1 induced the concentration of electron density between the reducing Lg n layers. The localization of high electron density on Doxorubicin in D-CNDs3 and D-CNDs4 suggested that it can react more easily in these two configurations than in the others.

Computational Methods
The Doxorubicin, CND and D-CND structures were drawn using Chem Sketch (ACD Lab, Toronto, CA, USA), and they were modelled using ArgusLab 4.0.1 [62,63] software without considering the discrete solvent medium. The calculations were performed using a personal computer with the Windows 10 operating system (built 19045) equipped with an Intel(R) Core(TM) i5-10210U CPU @ 1.60 GHz, 2112 Mhz, 4 cores, 8 logic gates and 16 GB of RAM installed. Conformational analysis was carried out through geometrical optimization using the PM3 semi-empirical quantum mechanical parameterization of the Hartree-Fock calculation method.
The geometry of Doxorubicin, CND and D-CND structures were considered optimized when the converged set point of 0.001 kcal/(Å mol) was reached using the Polak-Ribiere conjugate gradient algorithm for the optimization process [64] using up to 50,000 cycles. The interaction energies between the Doxorubicin and CNDs were calculated as follows: where Ei is the interaction energy, ED-CNDs is the total energy of the optimized Doxorubicin-CNDs structure, ED-CNDs is the total energy of the free optimized CNDs structure and ED is the total energy of the free optimized Doxorubicin. The optimized structures were also used for the numerical evaluation of the highest occupied molecular orbital (HOMO) and lower unoccupied molecular orbital (LUMO) energetic values. The electron density map showed a rich region between L g 1 and L g 2 in CNDs that was reduced in the D-CNDs, with the exception of D-CNDs2. In this case, the insertion of Doxorubicin between L h 1 and L g 1 induced the concentration of electron density between the reducing L g n layers. The localization of high electron density on Doxorubicin in D-CNDs3 and D-CNDs4 suggested that it can react more easily in these two configurations than in the others.

Computational Methods
The Doxorubicin, CND and D-CND structures were drawn using Chem Sketch (ACD Lab, Toronto, CA, USA), and they were modelled using ArgusLab 4.0.1 [62,63] software without considering the discrete solvent medium. The calculations were performed using a personal computer with the Windows 10 operating system (built 19045) equipped with an Intel(R) Core(TM) i5-10210U CPU @ 1.60 GHz, 2112 Mhz, 4 cores, 8 logic gates and 16 GB of RAM installed. Conformational analysis was carried out through geometrical optimization using the PM3 semi-empirical quantum mechanical parameterization of the Hartree-Fock calculation method.
The geometry of Doxorubicin, CND and D-CND structures were considered optimized when the converged set point of 0.001 kcal/(Åmol) was reached using the Polak-Ribiere conjugate gradient algorithm for the optimization process [64] using up to 50,000 cycles. The interaction energies between the Doxorubicin and CNDs were calculated as follows: where E i is the interaction energy, E D-CNDs is the total energy of the optimized Doxorubicin-CNDs structure, E D-CNDs is the total energy of the free optimized CNDs structure and E D is the total energy of the free optimized Doxorubicin.
The optimized structures were also used for the numerical evaluation of the highest occupied molecular orbital (HOMO) and lower unoccupied molecular orbital (LUMO) energetic values.
We also report the electron density map produced by the calculation ran using Argus-Lab 4.0.1 with a 40-slice grid.
The visualization of the molecules and geometrical parameters of the optimized structures were evaluated using Avogadro 1.2 software [65]. The layer distances are reported as the maximum distance, and the angular strain is defined as the dihedral angle formed by the layers.

The Interpretation of CDs Structures and Their Simulations: A Critical Discussion
The use of a computation model for the interpretation of CDs can be a powerful tool, but it should be used with great attention. The first limitation that the computational routes showed is related to the nature of what is simulated. As reported by Zhou et al. [54], the simple synthesis of CNDs led to the production of mixtures of great complexity, even using high-performance purification systems such as dialysis. The authors identified at least five well-characterized species inside the CNDs, with different distributions of functionalities. The authors proposed a data-driven interpretation model in which a layered structure was highly functionalized and variable. Nevertheless, the authors clearly stated that each fraction produced was not composed of a single specie but of a narrow distribution of very similar species. According to this capital finding, the computational approaches to CDs interpretation are intrinsically limited due to the mismatch between the simulated input and real materials. Nonetheless, there is a limited number of cases in which the gap between reality and the simulation is narrow enough to be neglected. GQDs can be investigated using the computation tool together with a complex set of electron microscopy and spectroscopic techniques, allowing for good agreement between the empirical data and computational outputs. On the contrary, PCDs are out of the range of possibility due to their random shapes. CNDs lay in between PCDs and GCDs, showing simple features to be integrated in a simulation such as a heptazine structure, and in more complex ones, distribution functions, as reported in Figure 5. We also report the electron density map produced by the calculation ran using Ar-gusLab 4.0.1 with a 40-slice grid.
The visualization of the molecules and geometrical parameters of the optimized structures were evaluated using Avogadro 1.2 software [65]. The layer distances are reported as the maximum distance, and the angular strain is defined as the dihedral angle formed by the layers.

The Interpretation of CDs Structures and Their Simulations: A Critical Discussion
The use of a computation model for the interpretation of CDs can be a powerful tool, but it should be used with great attention. The first limitation that the computational routes showed is related to the nature of what is simulated. As reported by Zhou et al. [54], the simple synthesis of CNDs led to the production of mixtures of great complexity, even using high-performance purification systems such as dialysis. The authors identified at least five well-characterized species inside the CNDs, with different distributions of functionalities. The authors proposed a data-driven interpretation model in which a layered structure was highly functionalized and variable. Nevertheless, the authors clearly stated that each fraction produced was not composed of a single specie but of a narrow distribution of very similar species. According to this capital finding, the computational approaches to CDs interpretation are intrinsically limited due to the mismatch between the simulated input and real materials. Nonetheless, there is a limited number of cases in which the gap between reality and the simulation is narrow enough to be neglected. GQDs can be investigated using the computation tool together with a complex set of electron microscopy and spectroscopic techniques, allowing for good agreement between the empirical data and computational outputs. On the contrary, PCDs are out of the range of possibility due to their random shapes. CNDs lay in between PCDs and GCDs, showing simple features to be integrated in a simulation such as a heptazine structure, and in more complex ones, distribution functions, as reported in Figure 5.  [54].
As clearly emerged, the Lh n layers have been defined, and their identification can be made using a combination of NMR and mass spectroscopy investigation. The definition of the Lg n layers is practically impossible and there is no solution to their unique identification. Furthermore, the interaction between Lh n and Lg n is only considered to be a weak interaction, but it is impossible to exclude the formation of a chemical bond between them, which could probably form as a consequence of the opening of epoxy functions. Therefore, the trade-off between representability and the reality of CND models should be evaluated considering how the model is able to predict the properties, but also, this approach could be misleading because the distribution of the particles could show the same properties as As clearly emerged, the L h n layers have been defined, and their identification can be made using a combination of NMR and mass spectroscopy investigation. The definition of the L g n layers is practically impossible and there is no solution to their unique identification. Furthermore, the interaction between L h n and L g n is only considered to be a weak interaction, but it is impossible to exclude the formation of a chemical bond between them, which could probably form as a consequence of the opening of epoxy functions. Therefore, the trade-off between representability and the reality of CND models should be evaluated considering how the model is able to predict the properties, but also, this approach could be misleading because the distribution of the particles could show the same properties as one of the particles. Nowadays, the computational study of CNDs could provide interesting insight into how these species interact with other molecules, but they are still far from being used as a predictive instrument with which to rationally design the synthesis and applications.
The customizable design of CNDs is of utmost relevance for all of those applications that require the fine tuning of the bot size and chemical functionalities [66].
Gao et al. [67] studied the use of the ab initio approach to CND hybrid materials coupled with CDs while assuming a layered structured for species. The authors showed the relevance of both the size and the shape of the layers in driving both of the relative positions and the alignment of bands. The authors proved that a considerably small enlargement size of the CNDs was able to enhance the visible light response of the species, forming a proper type-II van der Waals heterojunction between the CNDs and CDs [68]. The tuning of the band gap is also crucial for all of those catalytic applications, ranging from electrochemical to photochemical ones [69,70].
Li et al. [71] evaluated the effect of CNDs on the expression of protein kinases for the regulation of the cell signaling pathway. The authors showed the crucial role played by the phosphorylation of CNDs and how this is related to the edge functionalities.
The results discussed showed that the cutting edges applications of CNDs required a proper integration of empirical data with a solid interpretation that only a computational approach could provide. Nevertheless, some key issue remain unsolved together with opportunity as summarized in Table 2. Table 2. Summary of advantages and unresolved issue related to CNDs modelling for customdesigned applications.

Advantages Disadvantages
Comprehensive interpretations of the spatial arrangement of the CNDs structures. Evaluation of band gap. Description of electrochemical properties. Model the interaction between CNDs and biological species such as protein and nucleic acids.
Required a great analytic efforts to define the chemical functionalities. Limitation to GQDs and CNDs. Required long times. For several application is necessary only the knowledge of physiochemical properties.

Conclusions
The computational studies of CDs have made the first steps and already provided key information that have driven the research in new directions. Here, we describe a simple, representative and solid model of CNDs useful for evaluating their interaction with Doxorubicin, a widely used anticancer molecule. The modelled D-CNDs species showed various unique structures, suggesting that the intercalated species are more stable if the intercalation occurred between L g n or through adsorption directly onto L h n . Interestingly, π−π interactions seemed to play a greater role compared to hydrogen bond, but this result is quite sensitive to the model used. Nevertheless, this first computational study on threedimensional D-CND models enlightens the critical role of the modelling approach in order to fully understand the CDs' physiochemical behavior.
In the future, we aim to refine this model using a synthetic approach oriented towards the production of well-defined CNDs with reproducible and uniquely identifiable features. This will allow us to definitively prove the robustness of the approach in the prediction of the CDs' chemical, optical and geometrical properties.